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When simulating multicomponent mixtures via the Lattice Boltzmann Method, it is desirable to 
control the mutual diffusivity between species while maintaining the viscosity of the solution fixed. 
This goal is herein achieved by a modification of the multicomponent Bhatnagar-Gross-Krook (BGK) 
evolution equations by introducing two different timescales for mass and momentum diffusion. 
Diffusivity is thus controlled by an effective drag force acting between species. Numerical simulations 
confirm the accuracy of the method for neutral binary and charged ternary mixtures in bulk conditions. 
The simulation of a charged mixture in a charged slit channel show that the conductivity and electro- 
osmotic mobility exhibit a departure from the Helmholtz-Smoluchowski prediction at high diffusivity. 

PACS numbers: 


I. INTRODUCTION 

In recent years, multicomponent transport is enjoying growing interest due to burgeoning applications in biology, 
environmental science and energy management. One of the most ubiquitous transport phenomena regards electrolytic 
solutions and involves multiple processes, including fluid flow, mulfi-species diffusion and elecfrosfafic inferacfions. 
A better understanding of electrokinetic flows at nano and microscale level is paramount. At the nanoscale, for 
instance, it helps predicting mass and charge transport in biological ion channels. At the microscale, it guides 
the design of devices for biomolecular diagnostics and energy transfer systems, like micro fuel cells and batteries. 
Electrokinetic flow is also important for non-mechanical actuahng techniques, such as for pumping, mixing and 
separafing techniques ||l}j3l. 

Computer simulation is a direct approach to investigate microfluidic and nanofluidic systems and is exploited at 
best by those computational methods tailored to determine the concentration and current fields in generic confining 
geomefries. Among fhese, the Lattice Boltzmarm Method (LBM) is a reference fechnique 111151 and the reasons for ifs 
success are numerous: i) LBM is based on kinefic fheory and freafs the elementary interactions between particles via 
either a lumped model or top-down approach or by resolving the collisional kinetics at microscopic, bottom-up level, 
ii) the field dynamics is encoded by the single-particle distribution function and the governing equation is solved 
at second-order space /time accuracy by a compact numerical kernel, and iii) being the support a simple cartesian 
mesh, the method is amenable to efficient and parallel implementation, with complex boundary conditions handled in 
manageable and effective terms. In general, the evolution of the distribution frmction can be driven by several different 
interaction terms. Consequently, the hydrod 3 mamic fields can reproduce a vast niunber of physical conditions. In 
electrokinetics the mixture is composed by a neutral solvent and at least two charged species, with the Poisson 
equation solved numerically for electrostatics under appropriate boundary conditions. When ions are dissolved 
in water, they display a diffusion coefficient D ^ 10“^m^/s which is much smaller than the kinematic viscosity of 
the solution (at ambient temperature pure water has v = 10~^m^/s). In particular, charge transport displays an 
Ohmic contribution directly proportional to the diffusion coefficient, and a convective, electro-osmotic confribution 
inversely proportional to the kinematic viscosity. In order to accoimt quantitatively for the two contributions, it is 
crucial that the numerical method reproduces such disparity in the transport coefficients. A related problem emerges 
when modeling the different timescales involved in mass and momentum diffusion, where in a liquid the momentum 
field around a molecule diffuses much faster than the molecule itself. The Schmidf number, being fhe rafio befween 
kinemafic viscosify and fhe diffusion coefficienf (Sc = v/D), is as large as 10^ — 10^ in aqueous solufions. 

Perhaps the most popular implementation of fhe LBM is based on a reducfion, proposed about sixty years ago by 
Bhatnagar-Gross-Krook (BGK) jS), of the Boltzmann transport equation. The BGK scheme replaces the collision kernel, 
involving integrals of products of the phase space distribution fimchon, by a much simpler form, describing the 
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process of relaxation towards local equilibrium. Such an approximate treatment has certainly many advantages, both 
from the analytical point of view, because it allows to easily derive hydrodynamic equations and closed expressions for 
the transport coefficients, and from the computational point of view |7] . The BGK method, however, has an important 
limitation due to the fact that it is intrinsically built as a single-time relaxation process. As a consequence, the values of 
the mutual diffusion coefficient and kinematic viscosity are exactly equal. Several authors have proposed alternatives 
to lift such restriction by proposing, within the framework of a simplified kinefic model, failored inferacfions among 
particles of different species ISl- ITn . but only few of them have addressed the practical problem of fhe d 5 mamics of 
fluid mixfures, aparf from some nofable exceptions IIT2HT51 . 

Hamel proposed to control cross collisions by both an internal coupling force, proportional to the diffusion 
velocify, and an additional coupling ferm in fhe effecfive sfress fensor ll9l[T0l . Unforfunafely, fhe mufual diffusivify 
and the mixture kinematic viscosity cannot be independently controlled within a single cross-collision relaxation 
time. Subsequently, a multiple-relaxation-time (MRT) approach was proposed to independently control the mutual 
diffusivity and the mixture kinematic viscosity, the latter differing from the elementary mass averaged kinematic 
viscosity M- However, if a BGK-like equation for each species is assumed, the equations does not recover the 
single component dynamics. To overcome such drawback, Asinari proposed an alternative version within the MRT 
framework, bridging fhe MRT scheme wifh fhe Hamel model IITSlI . 

In fhe presenf paper, by using microscopically mofivafed physical argumenfs, we show how fo fake info accounf 
fhe occurrence of two different time scales associated with concentration and viscous momentum diffusion. Such 
a separate control was previously achieved by considering the fully microscopic approach, where the interactions 
between the particles were treated as hard-sphere collisions IITtIITsI . However, the fully microscopic approach can be 
exceedingly demanding in computational terms and unneeded in practical situations. Therefore, we propose a simpler 
version, which can be viewed as a reduction of the full method, by using an effective drag force exerted between 
fluid componenfs. Inspired by kinefic densify functional fheory lfT9l , we firsf infroduce a microscopic descripfion of 
fhe sysfem, in which a free paramefer defermines the diffusion coefficienf. We then extend the description to 
ternary charged mixtures, generalizing the self-consistent d 3 mamical method from neufral binary mixtures Il20l to 
electrolytic solutions. 

The paper is organized as follows: in secfion|^sfarfing from a Bolfzmann-like fransporf approach we infroduce a 
separation of fhe collision term into a BGK-type relaxation contribution plus a force term associated with the drag 
forces between differenf fluid componenfs which leads fo a modified BGK mulficomponenf equafion. Affer deriving 
fhe balance equations for fhe momenfa and fhe densifies, we show analytically thaf the diffusion and fhe viscosify 
coefficient display the required behavior. In section |III A we specialize the description to a neutral binary mixture, 
whereas in section III B|we consider a ternary charged mixture and characterize their transport coefficients. In section 
we corroborate the results by studying numerically the properties of the systems discussed above. In section[V|we 
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make some conclusive remarks and considerations. 


II. EFFECTIVE TREATMENT OF THE KINETIC EQUATION 


The exact time evolution of the one particle distribution function of the a species, /“(r,v, f), in a fluid mixture 
characterized by hard sphere interactions plus Coulomb forces can be formally written as 


|r(r,v,f)+v.Vr(r,v,f) + ?^ 


-r(r,v,f)=^n*/^(r,v,f) 


( 1 ) 


where jn'* are fhe particle masses, F“ an external force ferm and n'*^(r, v, t) describes fhe effecf of interaction between 
particles of type a and fi on the evolution of /“. When studying ternary charged mixtures, for instance, the index a 
can take the values , 0, -|-} depending on the particle charge. We assume that the interaction potential between two 
particles can be separated into a short-range strongly repulsive and a long range contribution (Ir) and correspondingly 
we make the approximation of splitting fhe collision ferm as: 

n'*^(r,v,f) = -in“/(r,v,f). (2) 

The repulsive part of , assimilable to a chard-sphere potential, could be modeled using the Revised Enskog theory 
(RET) of Ernst and van Beijeren ||2TI for hard-sphere mixtures which neglects velocity correlations for two particles 
about to collide, as in Boltzmann theory, but includes the configurational correlations resulting from the finite size 
of fhe particles via fhe pair correlafion function af confacf. However, even fhe solufion of RET equafions represenfs 
a formidable numerical problem, so fhat a simplified treatment is a valid alternative. In fact, in the present paper 
we shall not treat in detail the hard sphere collision term, as we did in previous papers IITtIITsII , and replace it by the 
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much simpler BGK relaxation term. The electrostatic contribution to the kinetic equation are treated in the mean field 
approximation, also known as Random Phase Approximation 1221 , so that to rewrite eq. 0 


+ v 


Vr(r,v,t) + 


F'*(r) d 




-r(r,v,t) = ^n“ 


(>-^v,t) + ^V</.(r).Ar(r,v,t) 


(3) 


where ip{r) is the self consistenf elecfric potential given by the solution of fhe Poisson equation 

V^ip(r) =- 47 r/B[n+(r) - n"(r)] , i/i = ei/^/keT 


(4) 


where Ig = /{4:TikgTe) is the Bjerrum length. The boundary conditions for a surface charge density E(r) are 

imposed as [Vi^(r)]j^ = —^^E(r), where the symbol _L indicates the gradient component orthogonal to the surface. 
A practical treatment of the collision operator was suggested by Dufty et al l23ll24l who, starting from the revised 
Enskog theory (RET) for hard spheres systems 1211 l25l . proposed to separate the contributions to stemming from 
fhe hydrod 3 mamic modes from fhe non-hydrodynamic ones. Such a goal is achieved by projecting fhe collision ferm 
onto fhe hydrodynamic subspace, sparmed by fhe frmcfions {l,v,v^}, and onto fhe complemenfary kinetic subspace: 


n-*/" = + {I- ■ 


(5) 


In the case of constanf temperafure, fhe projection of the collision operator onto the hydrod 5 mamic space is approxi¬ 
mated by the following formula (see l 20 l l 


Vi. j = 

' hydro^‘‘ 


^“(r,f) 


m^v 


Ta 


(v-u(r,f)) ■ <I>“(r,f) 


( 6 ) 


13 / 2 , 


being the local Maxwellian distribution and U 7 ^ = ^/k^T/nia 


with(/)'‘(r,v,t) = n“(r,t)[-^J-'-exp( 

^T,oc 

the thermal velocity of fhe a species. 

The terms ‘I>“(r, f) featuring in Eq. ||^ can be decomposed into different contributions stemming from different 
physical mechanism l26l : a term proportional to the gradient of the non ideal part of the chemical potential of species 
a, a drag force, a viscous force and a force proporfional to the gradient of the temperature. Eor the present goal, we 
shall only consider the drag force acting on the a species originating by the presence of the fi species for slow varying 
densities. This contribution is proportional to —?n“ 7 “^(u“ — u^) with 




» / m‘^m^ 2k bT 


where cr« is the molecular diameter of species a, = (na -|- cr ^)/2 and is the pair correlation at contact 

distance. In the case of molecules having approximately equal diameter cr and equal mass m, we can write 


8 9 IksT , , nP 


7 “^ = = <^drag'— 

3 V mTi° n ^ n 


(7) 


where 4*^® dimension of a frequency and will be used to time the diffusion coefficient in the model. ITereafter, 
we shall assume the following expression of fhe drag force : 


[u* (,, t) - uO (r, 01. 


m 


( 8 ) 


Concerning with the projection of 0^^ onto the non-hydrodynamics sub-space, Dufty and coworkers, who dealt with 
the one component case only, approximated it by a phenomenological single relaxation-time BGK prescription, which 
preserves the number of particles, the momentum and the kinetic energy and fulfills the physical symmetries and 
conservation laws of the fluid Il23l . Such a simple prescription, when extended to multicomponent fluids, has a serious 
drawback, as sfressed in IITSl . because if leads to the unphysical result that the diffusion coefficienf and fhe kinematic 
viscosity have the same value. In order to remedy such a situation we introduced the following approximation for the 
decay of fhe non hydrod 3 mamic modes: 


{I - V^yiroW 




(9) 
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The "orthogonalized" Maxwellian distribution is defined as 

... , K(r,f)-u(r,f)].?fi(v-u(r,f)) 

(p‘[{r,v,t) = (^'*(r,v,f)|H-2- 


■'Ta 


1 


+ ;r^[u'*(r,f) -u(r,f)][(u“(r,f) -u(r,f)] : '^^(v - u(r,f)) [, 


( 10 ) 


where are Hermite tensorial pol 5 momials of order k, such fhaf fhe kernel has a null projecfion onfo fhe {1, v,v^} 
subspace. The modified collision operafor (??) confribufes fo defermine fhe value of fhe shear viscosify, which furns 
ouf fo be a funcfion of ojyisc, a phenomenological collision frequency chosen fo reproduce fhe kinefic confribufion fo 
viscosify The final resulf is fhe following sef of Enskog-like equafions 


d 

dt 


r(r,v,f)+vVr(r,v,f) 


F*(r) 

m 


^r(r,v,f) 


F“'*«^(r,f) 


mv^ 


(v - u(r,f))</)“(r,v,f) 


3 

— V)/-(r).-r(r,v,f). 

m 3v 


( 11 ) 


Nofice fhaf in the case of a one-component fluid fhere is no difference between and (p“, since the velocities u® and u 
coincide. The above prescription fulfills fhe indifferentiabilify principle which sfafes fhaf when all physical properties 
of fhe species are idenfical, fhe fofal disfribufion f = musf obey fhe single species fransporf equation. 

The reason fo use fhe modified disfribufions in Eq. (T0| | rnsfead of Eq. ||^ is fo obfarn fhe correcf mufual diffusion 
and hydrod 5 mamic properfies sfarfing from Eq. 0- Erom fhe knowledge of fhe /®'s if is possible fo defermine nof 
only all fhe hydrod 5 mamic fields of inferesf, buf also fhe sfrucfure of fhe fluid af fhe molecular scale. The equafions 
consfifufing the building blocks of fhe classical elecfrokinefic approach can be derived from fhe sef of equafions <Tl||. 
In facf, fhe Poisson-Nernsf-Planck (PNP) and Navier-Sfokes (NS) equations Il27l are sfraighfforwardly recovered by 
faking fhe appropriafe velocify momenfs of fhe kinefic Enskog-Bolfzmann equafion. This derivation is performed in 
fhe nexf section in fhe limif of slowly varying fields. 


III. DERIVATION OF THE EQUATIONS OF ELECTROKINETICS FROM THE MICROSCOPIC APPROACH 
Sfarfing from fhe disfribufions functions /“ (r, v, f), we define fhe partial number densities, 

n“(r,f) = J dvf‘^{r,v,t), 


the mass densities p“(r, f), 
the species velocities u“(r, t), 

the barycentric velocity. 


^“(Tf) = ?naM“(r,f), 


u(r,f) = 


P(Tf) 


( 12 ) 


(13) 


(14) 


(15) 


with p(r, f) = p“ (r, f), and the kinetic contribution of fhe a componenf fo fhe pressure fensor 

= m‘^ J dv{vi-Ui){vj (16) 

We now rnfegrafe Eq. 0 w.r.f. fhe velocify and, assuming fhaf fhe disfribufion funcfions /“(r,v, f) go fo zero 
sufficiently fast, obtain the conservation law for fhe parficle number of each species 


|n“(r,f) = -V- (n“(r,f)(u'‘(r,f)-u(r,f))-V- (n'‘(r,f)u(r,f))= 


(17) 
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In order to recover the PNP equation, we consider the momentum balance for the species a which is obtained after 
multiplying by v and integrating w.r.t. v: 


— [n‘^{r,t)uJ{r,t)]+Vi(n‘^{r,t)uf{r,t)uJ{r,t) - n'*(r, f)(wf (r, f) - M,(r,!))(«“ (r, f) - Wy(r,f))): 

F“(r) F“- °(r.t) 

- -h - -M (r,f) + — 

jjja m"" ' ' ma 


- V, 


fr 

^n‘^(r,f)-^n''(r,f)V,,/;(r,f), 


(18) 


One can verify that the presence of fhe collision term (??) does not affect explicitly the hydrodynamic balance equations 
(??) and (??), thus allowing the decoupling of fhe concenfrafion diffusion from fhe momenfum diffusion which is fhe 
goal of fhe presenf work. If is convenienf fo infroduce fhe chemical pofenfial of fhe individual species, by fhe 
following equalify: 

n'‘(r,f)V,^“(r,f) = Vy7r“.(r, . (19) 


To derive fhe fofal momenfum equafion, we sum Eq. (??) over all componenfs and obfain fhe following expression: 


1 1 

dtUj{r:, t) + M,(r, f) V;My(r, t) + i^ij + - 


ez‘*n“(r,f)V^!p(r,f) 


1 

P 


Y n^{r,t)Ff{r)=0, 

oc=0,± 


( 20 ) 


where we have rnfroduced fhe fofal pressure fensor P,y(r, f) = YLa 0- Nofice fhaf fhe sum of fhe drag forces 
vanishes in (??) by fhe fhird principle of D 3 mamics. By neglecfrng shorf range rnferacfions, fhe pressure fensor is of 
kinefic nafure only and can be casf in fhe following form: 




\ 

dXi 3 dX]( ‘0 


( 21 ) 


where fhe diagonal parf is fhe ideal gas pressure of fhe mixture P,^ = n“ and t] is the dynamical shear viscosity. 

In the absence of elecfric and exfemal fields Eq. (??) is idenfical fo fhe one describing a one componenf fluid. Thus, 
by applying a sfandard Chapman-Enskog analysis, nof reporfed here for brevify Il28l , if is possible fo show fhaf fhe 
viscosify is relafed fo fhe relaxafion paramefer o^jj/sc/feafuring in Eq. | (Tl| |, as: 

(22) 

^visc OL 


a resulf applying fo mixfures of arbifrary number of componenfs. 


A. Diffusion and viscosity of a binary neutral mixture 


Let us specialize the treatment to two neutral species A and B, so that ip = 0, and the drag force achng on fhe A 
species due fo fhe B species is 


F^'‘^™^(r,f) 


n®(r,f) 




^drag 


assuming equal masses, negligible variafion of fhe densifies and negligible non linear ferms, Eq. (??) reads 
: (u^ (r, f) - u® (r, f)) + cv,rag (u^ (r, f) - u® (r, f)) = -1 { (r, f) - V/T® (r, t) - F^ (r) + F® (r)} . 

By neglecfing fhe relafive accelerafion, we obfain 


3 


u^(r,f)-u®(r,f) = - 


1 1 


m 


^drag 


{V;^^(r,f) - V/(r,f) -F^(r)+F®(r)}, 


that, in absence of external forces (F^ = F^ = 0), becomes 


u^(r,f)-u®(r,f) = - 


1 


^^drag 


(V/(r,f)-V/(r,f)). 


(23) 


(24) 


(25) 


(26) 
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Assuming also that the barycentric velocity vanishes, we rewrite Eq. (??) as 


h^(r,<) = V. 




In^+n^ mcOdrag 


(v?/^(r,t)-V/(r,t) 


(27) 


If the total density = uq is constant, so that = criQ and n® = (1 — c)mo are approximated by 

their ideal gas expressions, one finds 

|c(r,f)«^V2c(r,f) (28) 

^drag 

SO that the mutual diffusion coefficient is 



^drag 


(29) 


B. The charged mixture: conductivity and mass flow 


Let us turn back to the ternary charged mixture. Assuming that a steady current exists, we drop the non-linear 
terms in the velocities in the l.h.s. of Eq. (??), in absence of other external forces = 0, the approximated force 
balance, obtained from Eq. (??), reads 




(30) 


The r.h.s. of Eq. (??) represents the drag force exerted on the particles of t 5 q)e a = ±, in reason of their different drift 
velocities. In dilute solutions the charged components are expected to experience a large friction arising only from the 
solvent while a negligible friction from the oppositely charged species, so that we further approximate 


_ r^o 

^drag 


m 


l " -u°(Lf)) + "„(rY)^ (u^(r,f) -u^(r,f)) - -a;d™^[u^(r, f) - u(r, f)] (31) 

where we set n rP, that is, „± << „o 

. We also have that 

/±(r,t) = M±(r,f)u±(r,f) = M±(r,f)(u±(r,f) - u(r,f)) + M±u(r,f) « M±(r,f)u(r,f)--M±(r,f)F±'‘^™^(r,f). 

mCOdrag 

(32) 

Eliminating in Eq. (??), the ionic current in the stationary state is: 


- — - —M±(r,f)VM±(r,f)- - — 

mCVdrag mCVdrag 


r(h0 = - 


ez^n^ {r, t)Vip{r, t) -|-n*(r, f)u(r, f). 


(33) 


Eq. (??) is the phenomenological Planck-Nernst current, which is the sum of diffusive, migration and convective terms. 
The total electric charge density current is 


Je = - -M*(r,f)V|U*(r,f) -hUe/E-He^z^M*(r,f)u(r,f), 




'^drag 


(34) 


where the zero frequency electric conductivity Ug; is given by the Drude-Lorentz like formula 

aei = ^ ((z+)^M+ -h {z~fn-^ . (35) 

This expression shows that the conductivity is modulated by collisions with the solvent and decreases as the solvent 
becomes denser (cvdrag is an mcreasmg function of rf’) while increases with the number of charge carriers. Eor 
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completeness we write the macroscopic equation describing the electro-osmotic flow in the x-direction in a slit channel 
whose walls are normal to the z-axis, the so-called Stokes-Smoluchowski equation 


V xP CE; 


dz^ 


d^l 




dz^ 


= 0 


(36) 


and the Gouy-Chapman equation for the electric potential, which is based on the hypothesis that the charge density 
profile is in equilibrium in fhe direction orfhogonal fo fhe fluid flow. 


dz^ 


671 

ip(z) = — sinh 
e 


ksT ■ 


(37) 


In fhe nexf section we shall compare fhe analytic predictions of eqs. < [3^ and (??) wifh fhe numerical solutions of fhe 
kinetic model. 


IV. RESULTS 

We performed simulafions on binary neufral and femary charged sysfems based on fhe presenfed mefhod. In fhe 
discrefe LBM represenfafion, fhe sfandard procedure shows fhaf fhe kinetic equafion encoded by eq. (ITT) reads 1291 

/p(r + Cp,f-hl) = - cv^isc[fp{r,t) - (38) 

where a unif fime sfep is used and fhe vecfors {cp} are a sef of discrefe velocities used fo sample fhe disfribufion in 
special poinfs of velocify space. For fhe latter sef, we employ here fhe so-called D3Q19 velocify discretization l5l . 

The second ferm on fhe righf side of Eq ( [3^ describes fhe relaxafion towards fhe modified local equilibrium j^, 

while gp(r, f) represenfs the contribution due to the forces acting on species a. As shown in the Appendix, both terms 
arise from the Hermite expansion of the collisional kernel and are taken up to second order truncations of the Hermite 
series representation. The Chapman-Enskog analysis shows that the mutual diffusion coefficienf and fhe kinemafic 
viscosity are related to the input relaxation frequencies as D = — 1/2) and v = Vj{l/cOjjisc — 1/2), 

respectively In fhe simulafions, we always sef v = 1. 

We preliminarily fesfed fhe analytical predictions of Eq. (??) for fhe mufual diffusion coefficienf and of Eq. (??) for 
fhe electric conductivity in bulk conditions. In these tests we used x Ly x = 50 x 50 x 50 mesh points and 
applied periodic boundary conditions. To calculate the mutual diffusion coefficient of the binary mixture we prepared 
the uniform system by adding a small sinusoidal concentration profile, and monitored its exponential decay, having a 
wavelength dependent characteristic time 1/T{qz) = Dq^. The diffusion coefficienf was obfained a funcfion of 
and reporfed in Eig. The agreemenf wifh fhe fheorefical values is excellenf and fhe mefhod allows fo appreciably 
decrease fhe diffusion coefficienf. To compute fhe kinemafic viscosify we induced a sinusoidal shear modulation of 
fhe barycenfric velocify, and monifored ifs decay. We found fhaf upon changing fhe mufual diffusion coefficienf, fhe 
kinemafic viscosify remains unaltered, as expecfed. 

We next turn our attention to the charged ternary mixture. We first conducted a test to compute the electric 
conductivity in the bulk system subject to an uniform electric field. The density of the electrolytes is taken to be 
n^/n^ = 10“^. We measured the total electric current at varying diffusivity and at constant electric field. The 
data reported in Eig. show that the current is proportional to the diffusion coefficient D, as predicted by the 
Drude-Lorentz formula Eq. (??). 

We nexf analyze fhe fransporf behavior of fhe elecfrolytic solufion moving in a slif channel having charged walls, 
under fhe action of a imiform electric field parallel to the walls. The flow is directed along the x axis, the charged 
plates are aligned along the xy plane, with surface area Lx x Ly, and fhe charged plates are separated along fhe z axis 
by Lz. The fofal mesh size is x Ly x Lz = 40 x 10 x 50. The surface charge of fhe walls is TjLxLy/e = —0.04. 

Erom the macroscopic arguments used to derive the Hemholtz-Smoluchowski theory 1501 one expects that the 
mass current, 1,„ does not depend on the diffusion consfanf. The bulk densities of fhe charged species are fixed by 

imposing the Debye length, since = 1/^ATzlp{n'^ + wjj"), and we considered the three values Ap/Lz = 0.4, 0.2 
and 0.1. 

The total mass flow rate in the stationary state is defined: 

lm=m [ [n+{r)u+{r) + n~{r)u~(r) + n°{r)u°{r)]dz = m [ n{r)ux{r)dz 

Jo Jo 


(39) 
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By neglecting small variations of the total density with respect to the bulk value, the theoretical mass flow rate is 
given by [301 


1 


th 

m 


1 r kuU/l - 

V k^i tanh(A:Dkz/2). 


(40) 


where ki) = 1/ Ad is the the inverse of the Debye length. The quantity is an increasing function of Ad, explaining 
the results shown in Fig.j^ However, in the same plot we also notice a dependence of the mass flow rate on cOdrag- 
This behavior is not accounted for by the Helmholtz-Smoluchowski theory, which in fact shows no dependence on 
the inter-species diffusivity However, it should be borne in mind that one of the key assumptions of the theory is 
that the barycentric velocity locally equals the velocity of the neutral species, a condition that can be violated in the 
general case. 

To verify such assumption, we report in Fig.|^the differential velocity profiles of the three species and find that, 
when COdrag is large, the velocity of the neutral species is indistinguishable from the barycentric one. Vice versa, when 
the drag force is small the two velocities substantially depart from each other. 

The electro-osmotic contribution to the charge flow rale displays a similar dependence on cOdrag- To this purpose, 
we first consider the total electric current 


I = 



[n+(r)«+(r) 


— n (r)M^ (r)]dz. 


(41) 


and decompose it as the sum of two contributions stemming from Ohmic conduction and charge convection, 
I = lohm + leo, with the latter defined as 


Tn — 



-n (r)]ux(r)dz 


(42) 


and again is the barycentric velocity along the flow direction. Similarly, the electro-osmohc contribution to 
conductivity is obtained by dividing the electro-osmotic current by the applied tension. 

In the linear (Debye-Huckel) approximation, the convective contribution can be computed analytically and reads 


jth _ _ 

^eo 




sinh(EDTz)+A: dEz sinh(A:DTz/2) 


j/Ed tanh^(A:DTz/2) L 4cosh^(A:DTz/2) cosh{kr)Lz/2). 


(43) 


The behavior of the curves in Fig.|^is understood by noting that the number of charge carriers in the mixture is 
inversely proportional to A^ Il30l . Again, the electro-osmotic charge flow turns out to be an increasing function of Ad 
as shown in Fig. but depends on the mutual diffusivity, in particular by displaying a slight increase with D. The 
reason for such b^avior is the same as the one mentioned above. 

Having different velocihes between the solvent and the barycentric one at small cOdrag) that is at high diffusivity, 
gives rise to a larger barycentric velocity. We performed a theorehcal analysis of this occurrence at the level of 
linearized hydrodynamics starting from the kinetic equation (TT| |. The results show a fairly good agreement with the 
simulations data and predict that the phenomenon is only visible at small values of cOdrag- For the sake of conciseness, 
we will report the calculations elsewhere. 

In essence, the Lattice Boltzmann algorithm, derived from the kinetic model introduced in this work, is capable 
not only to accurately reproduce the behavior predicted by the continuum equations of Sec. Ill but also to exhibit 
non-trivial features of electrokinetic transport arising from the interplay of diffusion, convection and electrostatics. 


V. CONCLUSIONS 

To summarize, we modified the BGK d 5 mamics in order to control separately the timescales associated to con¬ 
centration and momentum diffusion. The presented method is derived from microscopic considerations involving 
the exact treatment of the collision term in the Boltzmann equation. We replaced the microscopic expression for the 
inter-species frictional force, which tends to equalize the species velocities, by a simpler effective term depending on 
a tunable frequency cOdrag- The derived equahons have similarities with the ones put forward by Luo and Girimaji 
M- However, the orthogonalization procedure described here represents a key ingredient to separately control 
concentration and momentum diffusion. 

We further implemented the equations in the context of the Lattice Boltzmann Method, and numerically verified 
that controlling independently the diffusivity and the kinematic viscosity for different physical models is feasible and 
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mirrors the expected dynamics in the continuum. We applied the method to the numerical calculation of the diffusion 
coefficient in neutral binary mixtures and of the electric currents of ternary charged mixtures under uniform bulk 
condifions, finding a perfecf agreemenf wifh fhe fheorefical prediction. 

Under non imiform conditions, such as those realized in a charged slit geometry, the electro-osmotic current exhibits 
a non-trivial behavior, since the convective contributions to the mass and charge currents display an interesting 
dependence on the diffusion coefficienf. Such occurrence, fo fhe besf of our knowledge, was nof previously reporfed 
in fhe liferature. 

We conclude with a remark concerning the relation between the present approach and the standard single relaxation, 
BGK description of mixtures. When the value of the tunable parameter is chosen to be equal to co^is^, on 
physical grounds one does not expect to observe differences between fhe fwo methods. Not only we have confirmed 
numerically such an occurrence, but it we also demonstrated that the kinetic equation (lT| | in practice reduces to the 
standard BGK single relaxation form in fhis limif. 
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Appendix A: Hermite expansion of the kinetic equation (TT) 

The standard procedure to construct the LBM algorithm starting from the kinetic equation is to expand each term 
of the equation in the Hermite polynomial basis. It is apparent that the coefficients of the expansion of the distribution 
function and the collisional terms are combination of the macroscopic fields Il29l . This property, together with the 
orthogonality of fhe basis set, not only guarantees that, given an expansion order K, the first K moments of the 
distribution are untouched by the truncation but also permits to control the error introduced by the evolution equation 
in the macroscopic fields II291I32I . Using fhe Einsfein convenfion on fhe repeafed indices, fhe one particle distribution 
function expands as 


/(r,v,f) =ro(v) Y, ~ro(v) Y =/^(r,v,f) 


k=0 

i = {h.-v4}; f/ = {1,2,3} ro(v) = 


K 

L 

k=0 ‘ 


{2nvj) 


,2^3/2 


exp(-^) ; ^f(v) = ^^Vfro(v) ; (Al) 


ro(v) 


where Vj represents the k-th derivative with respect to Vi^,Vi^, and <5!- is a sum of products of n 5's, each 6 having 
one subscripf from fhe sef i and one from fhe sef j Il33l . From fhe definifion, if follows fhaf 

fli(r,f) = [ f{r,v,t)'Hi{v)dv=< > 


(A2) 

In the following, we will often make use of the identities '^?y(v) = 'H? (v — u) + u, (u — v)y + uy(u — v), + u,uy and 
consider the Hermite expansion up to second order of fhe equafion 


wifh 


—/*(r,v,f)+v-V/“(r,v,f) = 


-pa,drag(j U 

g‘"{r,v,t) = -^-2-^ ■ {v -u{r,t))(pcc{r,v,t). 

mvt 


Both terms in the r.h.s. of Eq. (??) have a dependence of the type (v — so that 

= ^ / = 0 

= 4 /(v-u))^*dv= 

T) 'T' J / 'U'T'; 


T i 




■'T k 


v,vy - v^j-3ij){v - u)k(p‘^dv = 4 


■^T k 


ijk 


(A3) 


(A4) 


and 


= / (v - u)i(v - u)jCj)‘^dv + Ui / (v - u);^'*dv = VjS- 


s}jk = J (viVy - v^ij) (v - u)i,(p‘"dv = + Uy4,) 

The force ferm has fhe following discrefe represenfafion 


(A5) 


f) 


■ (cp - u(r,f))(/)4r,cp,f) - M“(r)ro(cp) i^p) + ^ 


T ij 


■^T ijk 


= M''(r)rop 
^ g“p(rH) 


prfrflg . „ I 

-^ + 4 (■ Cp)(u ■ Cp) - z;2F“'^™4r, t) ■ u 


(A6) 
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where Fop = ro(cp). The analogous expansion of (p‘^ reads 

= </>“(r,v,f){l + 4(u“-u)-?fi(v-u) + ^(u“-u)(u“-u):^2(^_y)| 

V ^7^ 'j' ^ 

= </)“(r,v,f) + 4(u“-u) :</)“(r,v,f)^i(v-u) + ^K-u) :</)“(r,v,f)^2(^_y) (^7) 

vt^ 2vj 

The expansion of fo second order is well-known (see for example Il29l l, so that we report here only the projections 
of — u), — u) to the same order. The integrals involved in the calculation of the projection of 

(p‘^V}{\ — u) , were given before (see Eq. (|A4}) while for the remaining ones, we obtain: 


= J ^“[(v - u),(v - u)y - v^Sij]dv = 0 

sfjk = J v/:^'*[(v - u)Kv - u)j - vlSij]dv = 0 

= j (v/V)t - - u)dv = V^iSiiSji: + Si^Sji) 


(A8) 


Putting all together, 
cp;,^ = n^Top{l + 


- - + - 

Vj 2vj - 


2 

2(u“ ■ Cp)(u ■ Cp) - 2z;p(u'‘ ■ u) - (u ■ Cp)^ -|- z;p|u|^ -h (^(u“ - u) ■ Cp^ - Vj\u‘^ - u|^ | 


(A9) 
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FIG. 1: Mutual diffusion coefficient for a neutral binary mixture in bulk conditions vs 0Jdrag/‘^visc- Solid line: theoretical curve, 
D/v = {^1 ^drag ~ ^1'^)I{^1 ^visc “1/2), Symbols: LBM data. 
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FIG. 2: Electrical conductance of a ternary charged mixture in bulk conditions vs the mutual diffusion coefficient. The conductance 
is divided by the reference value ug = e(w+ + n^)vjLyLz/ExLx- Solid line; Drude-Lorentz formula, Eq. (??), symbols: LBM data. 
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FIG. 3: Mass current vs Ap for the charged slit channel obtained with LBM and computed via Eq. (??). The mass current is divided 
by the reference value l„,o = /v. The curves are a guide fo fhe eye. 
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FIG. 4; Difference between the barycentric velocity and that of the neutral species for three diffusion coefficients as obtained via the 
LBM. Inset: velocity of the neutral species (solid line) and the barycentric one (symbols) corresponding to the two limiting cases 
^drag/<^visc = 0-35 (circles) and 5.25 (squares) and for Ap /L^ = 0.1. 
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FIG. 5: Electric conductivity for the charged slit channel as obtained via the LBM for Ap / Lz = 0.4 (circles), 0.2 (squares) and 0.1 
(diamonds). The linear dependence of the electric current on diffusivity reflects the dominance of the Ohmic contribution, which is 
about 100 times larger than the convective counterpart. 
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FIG. 6: Electro-osmotic contribution to the electric conductivity (see Eq. (??) and text for details) for fhe charged slit channel for 
Ajj/Lz = 0.4 (circles), 0.2 (squares) and 0.1 (diamonds). 










